Defect Annihilation and Proliferation in Active Nematics 



Luca Giomi, 1 Mark J. Bowick, 2 ' 3 Xu Ma, 2 and M. Cristina Marchetti 2 ' 3 

1 SISSA, International School for Advanced Studies, Via Bonomea 265, 34136 Trieste, Italy 
2 Physics Department, Syracuse University, Syracuse NY 13244, USA 
3 Syracuse Biomaterials Institute, Syracuse University, Syracuse NY 13244, USA 

Liquid crystals inevitably possess topological defect excitations generated through boundary con- 
ditions, through applied fields, or in quenches to the ordered phase. In equilibrium, pairs of defects 
coarsen and annihilate as the uniform ground state is approached. Here we show that defects in 
active liquid crystals exhibit profoundly different behavior, depending on the degree of activity and 
its contractile or extensile character. While contractile systems enhance the annihilation dynamics 
of passive systems, extensile systems act to drive defects apart so that they swarm around in the 
manner of topologically well-characterized self-propelled particles. We develop a simple analytical 
model for the defect dynamics which reproduces the key features of both the numerical solutions 
and recent experiments on microtubule-kinesin assemblies. 



Active liquid crystals are nonequilibrium fluids com- 
posed of internally driven elongated units. The key fea- 
ture that distinguishes them from their well-studied pas- 
sive counterparts is that they are maintained out of equi- 
librium not by an external force applied at the system's 
boundary, such as an imposed shear, but by an energy 
input on each individual unit. Examples include mix- 
tures of cytoskeletal filaments and associated motor pro- 
teins, bacterial suspensions, the cell cytoskeleton, and 
even nonliving analogues, such as monolayers of vibrated 
granular rods [1 j. The internal drive that characterizes 
active liquid crystals dramatically changes the system's 
dynamics and yields novel effects arising from the inter- 
play of orientational order and flow, such as spontaneous 
laminar flow [2HI], large density fluctuations [5-7 , un- 
usual rheological properties [8 10], excitability [TTJ [12] 
and low Reynolds number turbulence [T2j [13] . 

Ordered liquid crystalline phases of active matter, like 
their equilibrium counterparts, exhibit distinctive inho- 
mogeneous configurations known as topological defects. 
In equilibrium, defect configurations may be generated 
through boundary conditions, externally applied fields, 
or rapid quenches to the ordered state. When the con- 
straints are removed or the system is given time to equi- 
librate, the defects ultimately annihilate [14 . Experi- 
ments have shown that in active systems, in contrast, 
defect configurations can occur spontaneously in bulk 
and be continuously regenerated by the local energy in- 
put p~5j[16]. The nature of the topological defects is, of 
course, intimately related to the symmetry of the sys- 
tem, which can be either polar (like in ferromagnets) or 
nematic. While the nature of the charge +1 defects that 
occur in polar active systems has been studied for some 
time [T7lf2Q] , the properties of defects in apolar or ne- 
matic active media are still largely unexplored. In these 
systems the defects are charge ±1/2 disclinations [21]. 
Such defects have been identified in monolayers of vi- 
brated granular rods [7 and also in active nematic gels 
assembled in vitro from microtubules and kinesins. In the 
latter case the defects were shown to drive spontaneous 



flows in bulk [TBI l22]- When confined at an oil- water 
interface, furthermore, the gel forms a two-dimensional 
active nematic film, with self-sustained flows resembling 
cytoplasmic streaming and the continuous creation and 
annihilation of defect pairs [T6] . 

In this Letter, we examine the effect of activity on the 
dynamics of disclinations in a two-dimensional nematic 
liquid crystalline film [23]. Hydrodynamics plays an im- 
portant role in controlling the dynamics of defects in liq- 
uid crystals. As the defect moves, the coupling between 
the orientational order parameter and the flow velocity of 
the fluid yields what is usually called the backflow which 
significantly modifies defect dynamics [24H28] , Here we 
show that active stresses dramatically affect the defect 
dynamics by altering the backflow in such a way as to 
speed up, slow down, or even suppress pair annihila- 
tion, according to the extent of activity and the typi- 
cal time scale of orientational relaxation of the nematic 
phase. Moreover, when the latter is very large compared 
to the time scale associated with activity, relaxation is 
overwhelmed entirely, leading to defect proliferation. 

The hydrodynamic equations of active nematic liquid 
crystals can be obtained from that of passive nematics by 
the addition of nonequilibrium stresses and currents due 
to activity [TJ [11] [12] . These equations are formulated 
in terms of a concentration c, a flow velocity v and the 
nematic tensor order parameter Qij = S [nirij — |<%), 
with n the director field. The alignment tensor Qij is 
traceless and symmetric, and, hence, has only two inde- 
pendent components in two dimensions. Considering for 
simplicity the case of an incompressible fluid of constant 
density p, where V • v = 0, the equations are given by 

Dc 

£^ = di [ D ijdjC + a^djQij] , (la) 

P 15t = T1 ^ 2vi ~ diP + dj(Tij ' ( lb ) 
j^r = ^Suij + Qik^kj — UikQkj + 7 -1 ^ij •> (1c) 
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FIG. 1: (color online) Snapshots of a disclination pair shortly 
after the beginning of relaxation. (Top) Director field (black 
lines) superimposed on a heat map of the nematic order pa- 
rameter and (bottom) flow field (arrows) superimposed on a 
heat map of the concentration for an extensile system with 
a = —0.2 (a),(c) and a contractile system with a = 0.2 

(b) ,(d). In the top images, the color denotes the magnitude 
of the nematic order parameter S relative to its equilibrium 
value So = y/l — c* Jcq. In the bottom images, the color de- 
notes the magnitude of the concentration c relative to the 
average value cq. Depending on the sign of a, the backflow 
tends to speed up (a > 0) or slow down (a < 0) the annihi- 
lation process by increasing or decreasing the velocity of the 
+1/2 disclination. For a negative and sufficiently large in 
magnitude, the +1/2 defect reverses its direction of motion 

(c) and escapes annihilation. 



where §- t = 
Dij = Do 5, 



dt + v • V indicates the material derivative, 
%3 + DiQij is the anisotropic diffusion tensor, 
r] is the viscosity, p is the pressure, and A is the nematic 
alignment parameter. Here Uij = (diVj + djVi)/2 and 
Uij = (diVj — djVi)/2 are the symmetrized rate of strain 
tensor and the vorticity, respectively. The molecular field 
Hij embodies the relaxational dynamics of the nematic 
phase (with 7 a rotational viscosity) and can be obtained 
from the variation of the Landau-De Gennes free energy 
of a two-dimensional nematic [21], = —5F/5Qij, with 

F/K = JdA [\(c-c*)trQ 2 + \c(trQ 2 ) 2 + \\VQ\ 2 } , 

(2) 

where K is an elastic constant with dimensions of en- 
ergy, trQ 2 = S 2 /2 and c* is the critical concentration for 
the isotropic- nematic transition, so that, at equilibrium, 



S — y/l — c*/c. Finally, the stress tensor cr^- = cr\ - + erf ■ 
is the sum of the elastic stress due to nematic elasticity, 



QikH kj - H ik Qkj, where for simplic- 



a lj ~ —^SHij -f i^ik-tikj — riik^kj, 
ity we have neglected the Eriksen stress, and an active 
contribution, afj = ot2(?Qi^ which describes contrac- 
tile or extensile stresses exerted by the active particles 
in the direction of the director field. In addition, activ- 



ity yields a curvature-induced current j 3 



-aic 2 V • Q 



in Eq. (|TJl) that drives units from regions populated by 
fast-moving particles to regions of slow- moving particles. 
The c 2 dependence of the active stress and current is 
appropriate for systems where activity arises from pair 
interactions among the filaments via cross-linking motor 
proteins. The sign of 0L2 depends on whether the active 
particles generate contractile or extensile stresses, with 
0L2 > for the contractile case and <%2 < for extensile 
systems, while we assume a± > 0. 

To study the dynamics of defects, we consider a pair 
of opposite-sign half-integer disclinations separated by 
a distance where x± is the x coordi- 

nate of the ±1/2 disclination, respectively, as shown in 
Figs. [TJa) andjjjb). When backflow is neglected, the 
pair dynamics is purely relaxational and is controlled 
by the balance of the attractive force between defects 
^pair = -V£;p a i r , with £ pair - K\og{x/a) the energy of 
a defect pair (with a the core radius), and an effective 
frictional force Ff ric = fix, with \i ~ 7 a friction coef- 
ficient. Thus fix = K/x and the distance between the 
annihilating defects decreases according to a square-root 
law, x(t) (x y/t a — t, with t a the annihilation time. More 
precise calculations have shown that the effective fric- 
tion is itself a function of the defect separation [29l [30] , 
fi = fio log(x/a), although this does not imply substantial 
changes in the overall picture. This simple model predicts 
that the defect and antidefect approach each other along 
symmetric trajectories. 

We have integrated numerically Eqs. ([T]) for an ini- 
tial configuration of uniform concentration and zero flow 
velocity, with two disclinations of charge ±1/2 located 
on the x axis of a square L x L box at initial positions 
#±(0) = (=LL/4, 0). The integration is performed using 
the finite differences scheme described in Ref. [TTJ [12] . 
To render Eqs. ([!]) dimensionless, we normalize distance 
by the approximate length of the active rods £ = 1/ \f& 
stress by the elastic stress of the nematic phase a = Ki~ 2 
and time by r = r]£ 2 / K representing the ratio between 
viscous and elastic stress. In these dimensionless units, 
for simplicity, we let (%2 = a and take a\ = |qj2 |/2. 
Periodic boundary conditions are assumed, and the de- 
fects are allowed to evolve until they annihilate. Figure 
[T] shows a snapshot of the order parameter and flow field 
shortly after the beginning of the relaxation for both a 
contractile and extensile system, with a = ±0.2 in the 
units defined above (see also the supplementary movie 
SI [31]). 

In passive nematic liquid crystals (i.e., a = 0) it is 
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FIG. 2: (color online) Defect pair production in an active sus- 
pension of microtubules and kinesin (top) and the same phe- 
nomenon observed in our numerical simulation of an extensile 
nematic fluid with 7 = 100 and a — —0.5. The experimental 
picture is reprinted by permission from T. Sanchez et al., Na- 
ture (London) 491, 431 (2012). Copyright 2012, Macmillan. 
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well known that the dynamics of defects is greatly mod- 
ified by the so-called backflow, that is, the flow induced 
by reorientation of the nematic order parameter through 
the elastic stresses a\- in the Navier-Stokes equation. In 
particular, when backflow is neglected, the defect and an- 
tidefect are predicted to move at the same velocity toward 
each other until annihilation. Backflow tends to speed up 
the +1/2 defect and to slow down the —1/2 defect, yield- 
ing asymmetric trajectories [25]. In active liquid crystals, 
the active stress in the Navier-Stokes equation provides 
a new source for flow associated with inhomogeneities 
in the order parameter, as demonstrated first in a one- 
dimensional thin film geometry where activity drives a 
transition to a spontaneously flowing state [2] . This new 
active backflow can greatly exceed the curvature-driven 
backflow present in passive systems. Furthermore, the di- 
rection of the active backflow is controlled by the sign of 
the activity parameter a and, for a given director config- 
uration, has opposite directions in contractile and exten- 
sile systems. Backflow arising from active stresses drives 
the +1/2 defect to move in the direction of its "tail" in 
contractile systems (a > 0) and in the direction of its 
"head" in extensile systems (a < 0), where the termi- 
nology arises from the cometlike shape of +1/2 defects. 
In contrast, due to symmetry considerations, the active 
backflow vanishes at the core of a —1/2 defect which thus 
remains stationary under the action of active stresses. 
We note that active curvature currents in the concen- 
tration equation controlled by ql\ have a similar effect, 
as first noted by Narayan, Ramaswamy, and Menon and 
collaborators in a system of vibrated granular rods 0. 
Such active curvature currents control dynamics in sys- 
tems with no momentum conservation but are very small 
here, where the concentration variations remain small, as 
seen from Figs. [TJc) and[TJd), and flow controls the dy- 
namics. 

In contractile systems active backflow yields a net 
speed-up of the +1/2 defect towards its antidefect for 
the annihilation shown in Fig. [ljb) . In extensile sys- 
tems, with a < 0, backflow drives the +1/2 defect to 



FIG. 3: (color online) Defect trajectories and annihilation 
times obtained from a numerical integration of Eqs. |l]) for 
various 7 and a values, (a) Defect trajectories for 7 = 5 
and various a values (indicated in the plot). The upper (red) 
and lower (blue) curves correspond to the positive and neg- 
ative disclination, respectively. The defects annihilate where 
the two curves merge, (b) The same plot for 7 = 10. Slow- 
ing down the relaxational dynamics of the nematic phase in- 
creases the annihilation time and for a — —0.2 reverses the 
direction of motion of the +1/2 disclination. (c) Defect sepa- 
ration as a function of time for a — 0.2 and various 7 values, 
(d) Annihilation time normalized by the corresponding anni- 
hilation time obtained at a — (i.e., t®). The line is a fit to 
the model described in the text. 

move towards its head, away from its —1/2 partner in 
the configuration of Fig. [ljb) , acting like an effectively 
repulsive interaction. This somewhat counterintuitive ef- 
fect has been observed in experiments with extensile mi- 
crotubules and kinesin assemblies [16] and can be under- 
stood on the basis of the hydrodynamic approach em- 
bodied in Eqs. ([!]). In Fig. [2] we have reproduced from 
Ref. [16 a sequence of snapshots showing a pair of ±1/2 
disclinations moving apart from each other together with 
the same behavior observed in our simulations. 

To quantify the dynamics we have reconstructed the 
trajectories of the defects by tracking the drop in the 
magnitude of the order parameter. The trajectories are 
shown in Figs. |3ja) and[3|b), where red lines in the up- 
per portion of the plots represent the trajectory of the 
+1/2 disclination, while the blue lines in the lower por- 
tion of the plot are the trajectories of the —1/2 defect. 
The tracks end when the cores of the two defects merge. 
For small activity and small values of the rotational fric- 
tion 7, the trajectories resemble those obtained in Ref. 
[25] for passive systems. At large values of activity, how- 
ever, the asymmetry in defect dynamics becomes more 
pronounced, and when the activity dominates over orien- 
tational relaxation, the +1/2 disclination moves indepen- 
dently along its symmetry axis with a velocity v ex — a x, 
whose direction is dictated by the sign of a. This be- 
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havior is clearly visible in Fig. [3^c), showing the defect 
separation x(t) as a function of time. For 7 sufficiently 
large, the trajectories are characterized by two regimes. 
For large separation the dynamics is dominated by the 
active backflow, and thus x(t) oc —a and x(t) oc —at. 
Once the defects are about to annihilate, the attractive 
force F pair oc 1/x takes over, and the defects behave as 
in the passive case with x(t) oc y/t a — t. 

Building on these results, we now propose a phe- 
nomenological one-dimensional model that captures 
qualitatively the dynamics of pair annihilation in active 
nematics. By neglecting for simplicity the position de- 
pendence of the friction, which we assume constant, the 
dynamics of a pair of disclinations initially at a distance 
xq along the x axis is governed by the equations 



ja[x ± -v b (x±)] = =F- 



K 



x + 



(3) 



where v^{x±) is the backflow field at the position x± of 
the ±1/2 defect, given by v\>(x) = v+(x — x+) + V-(x — 
£_), with v±(x) the flow field due to an isolated ±1/2 
defect. We retain only the active contribution to the 
backflow and replace the flow profiles by their constant 
values at the core of the defect, with Vb(x+) = v a oc —a 
and Vb(x-) = 0. Note that v a < for contractile systems 
and v a > for extensile ones. This yields the following 
simple equation for the pair separation x = v a - 
with k = K/fjL. This equation explicitly captures the two 
regimes shown in Fig. |3jc) and described earlier. The 
solution takes the form 



x(t) = xq ± v a t — 2(n/v a ) In 



x(t) — 2tz/v a 



xq — 2k/v q 



(4) 



The pair annihilation time t a is determined by x(t a ) =0 
and is given by t a = -x /v a -(2K,/vl) In [1 - (x v a /2hi)]. 
For passive systems (a = 0) this reduces to t° a = x^/Ak. 
In contractile systems, activity speeds up pair annihila- 
tion, while it slows it down in extensile systems. Our 
simple model predicts that the annihilation time, nor- 
malized to its value in passive systems, t a /t®, depends 
only on v a xo/2K ~ cry. Figure |3jd) shows a fit of the 
annihilation times extracted from the numerics to this 
simple formula. The model qualitatively captures the 
numerical behavior. 

While the effect of activity on the precollisional dynam- 
ics of a disclination pair can be accounted for relatively 
simply in terms of active backflow, the postcollisional be- 
havior is dramatically affected by activity [32] . Figure [4] 
shows the evolution of the system after annihilation of 
the initial defect pair (see also the supplementary movie 
S2 [31]). The frame in Fig. |4ja) shows the initial configu- 
ration of the two defects, while FigjIJb) shows the config- 
uration just after pair annihilation. The other frames dis- 
play the evolution in time (with time increasing from left 
to right and top to bottom). Immediately after collision, 




FIG. 4: (color online) Schlieren texture highlighting the post- 
collisional dynamics of a ±1/2 pair for 7 10 3 and a 0.2. 
(a) indicates the initial configuration of the defects, and (b) 
shows the system immediately after defect annihilation. 



the system develops two density or flow bands reminis- 
cent of those observed in the absence of defects [12 . The 
bands, however, are unstable and quickly start deform- 
ing while new defect pairs "pinch off". The dynamics 
quickly becomes chaotic, with frequent defect formation 
and annihilation events in the background of an overall 
proliferation of defects. The passage of defects through 
a region of space lowers the local nematic order param- 
eter in that region. At large friction 7, the slow relax- 
ation prevents the restoration of the order parameter to 
its initial value, leading to a progressive reduction of the 
average order parameter in time. Complex textures in 
active nematics were also reported in Ref. [10 , although 
those structures are not easily decomposed in terms of 
disclinations. More work is needed to fully explore this 
rich and complex dynamics and formulate a quantitative 
classification of the behavior of defects in active liquid 
crystals. 
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